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Abstract. - Matsubara Green's functions for interacting bosons are expressed as classical sta- 
tistical averages corresponding to a linear imaginary-time stochastic differential equation. This 
makes direct numerical simulations applicable to the study of equilibrium quantum properties 
of bosons in the non-perturbative regime. To verify our results we discuss an oscillator with 
quartic anharmonicity as a prototype model for an interacting Bose gas. An analytic expression 
for the characteristic function in a thermal state is derived and a Higgs-type phase transition 
discussed, which occurs when the oscillator frequency becomes negative. 
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Introduction. - A commonly used method to study thermal properties of interacting many- 
body systems is the Matsubara Green's function technique |l], g, |J. In this paper, we present 
a method of evaluating Matsubara Green's functions for bosons, based on constructively 
characterising Feynman paths as solutions to an Ito stochastic differential equation (SDE). 
(For stochastic calculus and SDEs see M; a brief discussion is given below.) Numerically 
simulating this SDE then allows one to calculate physical quantities using their expressions as 
path integrals [bl (this constitutes a constructive rather than a Monte-Carlo method H). Since 
path integrals correspond to an in principle exact summation of the Matsubara diagram series 
as a whole (as opposed to approximate partial summations underlying conventional Green's 
function approaches P, H]), the proposed technique is applicable beyond the perturbation 
region and is only limited by numerical errors. It could be advantageous, e.g., for analysing 
the behaviour of an interacting Bose gas near the critical temperature of condensation M] . 

We develop a version of the Matsubara technique where a diagram series is averaged over the 
thermal P-distribution for free bosons M . For the grand canonical ensemble this distribution is 
positive and hence can be interpreted as a probability density. Then we see that summing the 
series is equivalent to solving a certain SDE. For bosons with quartic interaction, this SDE is a 
linear imaginary-time Schrodinger equation with multiplicative noise. As a result, we express 
normally ordered thermal averages of the bosonic field as classical statistical averages. We 
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verify this result, including the choice of stochastic calculus (Ito), by considering a prototype 
model of interacting Bose fields — a quantized oscillator with quartic anharmonicity A simple 
analytic expression for the characteristic function is derived and a Higgs-type phase transition 
discussed, which occurs when the frequency of the oscillator becomes negative. 

General theory. - Consider a system of interacting Bose fields in an external potential 
V = V(f), described by the Schrodinger field operator %[> = ip(r), with the Hamiltonian 
H = H + Hi n t, where 



H Q = J d 3 rft f-l-V 2 + Vji), H 



iat = # int (^,V0 = ~ /d 3 f^V 2 , (1) 



K = Aita/m and a is the s-wave scattering length. We use units where % = 1. 

The aim of the paper is to calculate thermal averages of normally ordered products of ijj 
and ffl over the grand-canonical density matrix, p = c( n ~ H+lJ - N ^ . Here /3 = 1/fcsT, T is the 
temperature, p is the chemical potential, and Z — exp{— f2/3} is the partition function. To 
this end, we consider the normally-ordered characteristic functional, 

X (£, D = (exp (tft) exp (-C4) ) , (2) 

where ^ijv = J d 3 r^(r)ip^(f) (say) and £ = £(r ) is an arbitrary complex function. 

To relate averages over p to those over the density matrix for the free system, po — 
e {n o -H o +nN)0 ^ we introduce the Matsubara interaction picture, 

^ T (r,T) = e ( H o-^)r^-{H -„N)r^ (3) 

} T (f,r) = e<- Ho -^ T ^e-( Ho -^ T = ^ T (r,- T ), (4) 

and the thermal "scattering matrix" , 



S = T r exp 



-H$ t ,$t) , (5) 



1-0 



where TL(ip T , ipr) = J_g dr7?i n t('0T: V't), and T T is the "time" -ordering operator with respect 
to r. (Note that we consider Matsubara evolution from r = —(3 to r = 0, not from to j3 
as usual: r = is the latest, not the earliest, "time" ; some of our relations therefore differ 
slightly from those in the literature.) Using p = Spo/(S)o, where (. . .)o denotes the "free" 
averaging over p , we find with £(£,£*) = x(£>D (^)o- 

X& O = (TVexp [-H(i> T , 4> T ) + e^r(+0) - r^r(-O)] ) Q • (6) 

Matsubara diagram series W, 0, 0] are obtained by expanding the exponent in (o) in a power 
series, changing the order of the summation and averaging, and then factorising the many- 
operator averages into products of free thermal Green's functions 11 ffl . Physical quantities are 
then obtained isolating certain "leading" classes of diagrams which may be summed yielding 
Dyson equations. This inevitably involves approximations which fail if the interaction is not 
small or in the vicinity of a phase transition. In order to have an approach which is in principle 
exact, we here proceed in a different way. We first apply Wick's theorem proper so as to bring 
the operator expression in (|6|) to normal order, and then employ the relation, 
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Here, F is an arbitrary functional of %j)(r ) and ^(f), fD 2 ifto is a functional integration and 
P(ipo) is the Glauber P-distribution for the non-interacting bosons in a thermal state, which is 
positive S. Wick's theorem can be given a compact form using functional derivatives Jul [ll| 



T r F('ip T ,'4> 1 



e A F( 



'4>—nl> T ;<t>— t ''l>T 



where 



A = 



dTidr 2 // d rid r 2 Z>r(ri, r 2 ; n — r 2 )- 



-/3 J J $(f>(ri,Ti)5(f)(r 2 ,T2) 

(r , r) and 6(r , r) are arbitrary independent c-number functions, and 



(8) 



(9) 



D T (fi,f 2 ;Ti -r 2 ) 



T T ?/i T (ri,Ti)^ T (f 2 ,r 2 ) 0) = 6»(ri - r 2 ) ^t(^*i, n), ^ T (r 2 , r 2 ) (10) 



is the Matsubara contraction. Note that Eq.(|h does not refer to any state vector and that 
Dt in (nO) is a vacuum and not a thermal average: as opposed to the usual Matsubara 
diagram series the thermal properties of the non-interacting system are contained in the initial 
distribution P(ipo) and not in the propagator. It is straightforward to prove that at the same 
time Dt is the retarded Green's function of the imaginary-time free Schrodinger equation: 



d_ 1_ 

v <9ti 2m 

Applying (||) to (||), and using (Q), we find 

x(£,t) ="*, * = e A exp[-H(0,0) +£0(0) - £*0(O) 



V? + V(n )-fx) D T {r u r 2 ; n - r 2 ) = 5(n - r 2 )5 3 (fi - r 2 ). (11) 



(12) 



Here ^o = V'o^j 1 ") and V'o — V'o^j 1 ") — V'cK^i - T ) are the amplitudes of the non- 
interacting fields in the coherent state |V>o) : iprir, r) |^o) = i>o(r, r) |^>o) an( i (V'ol ^"tC^i T ) = 
(V'ol V , o(' ? ) T )- F° r T = the Matsubara interaction picture coincides with the Schrodinger 
picture hence ^(^ , 0) = -00 if) and ^o(r , 0) = ipo{f). 

Our goal is now to express (O) as a classical statistical average. Since P(ipo) > 0, this 



interpretation applies to the averaging over -0o (denoted by the upper bar in (12)). In the 



quantity $ one can replace £0(0) by £i/>g, since r = is the largest "time" and the exponential 
operator does not act on 6(0) . Following the Stratonovich-Hubbard transformation used in 
path-integral approaches |l2j, we introduce a Gaussian stochastic variable r](r,T), such that 
T)(F,T)r)(f ,r') = «(5 3 (f — r )6(t — r'), and write, 

-H(4>, 



exp 



= exp 



\i(fa4j , 



where &r\& = J_ g dr J d 3 r 6{f , r) ^(f , r) 0(r , r). We then show that 



e A exp 



i6r)S - £*0(O) 



\tp—ti/; ;4>—tipo 



exp 



#0^ - CV'(O) 



where tp = tp(r, r) is a solution of the Ito equation, 



o -I 

— V 2 + V - \x ) 0(f , r) = w/(r , r)V»(f , r), 



with the initial condition given by 

^(f,-/3)=Vo(r,-/3). 



(13) 

(14) 

(15) 
(16) 



EUROPHYSICS LETTERS 



Indeed, expand all exponents on the LHS of ( pL4| ) in power series and perform the functional 
differentiations. The emerging series is similar to a diagram series for a problem with a 
linear perturbation. In particular, all connected diagrams are linear chains, of the structure 



either iPoDtitiDt ■ ■ ■ iTjif^o, or — £* DtitiDt ■ ■ ■ irjipo- By virtue of Meyer's first theorem |12 
this immediately results in the RHS of (UJ), with ijj obeying the Dyson equation, 



ip(f, t) = ipo(r,T) + I d r' / dr' Dt{t,t' ,r — r')ir](r' ,T')ip(f ,t'), (17) 

which in turn is equivalent to (fla). 

For r\ = 0, Eq.(|l5|) describes the evolution of f/'Oi while changing r — > — r yields the equation 
for tpQ. Hence -^ fd 3 fipo(f, T)ip(r, t) — i fd 3 ripo(r, T)rj(r, T)ip(f, r) and we find that 

o r 

3„- 



iiporjij} = i / dr \dri})o(f,T)ri{r,T)ip{r,T)= \d ripQ(f) [ip(r, 0) — tpo(f)] . 



Finally, 



x(£, D = m, D/x(o, 0), x(£, D = e A , (18) 

A= / d¥ [^W-r(W.0) + ^(W,l))- IV'oWI 2 ] . (19) 



The double upper bar in (JTq) denotes an averaging over both ipoi"?) and r](r,T); i.e., the 
statistics here is that of the random trajectories (= Feynman paths) tj)(r,T), which are 
solutions to the SDE (|l5[), with the random initial condition distributed with probability 
P(ipo). The interaction between the bosons is contained in the noise r\ while the thermal 
properties of the system in P(ipo)- Relation ( |lS| ) is the main result of this paper. 

It is important that its derivation implies a causal regularisation of the function Dt- 
Namely, Dt (which is a generalised function) is replaced by a certain number of times contin- 
uously differentiable function, yet preserving the causality condition DT(r,r , r) = 0, r < 0. 
Consider, for example, relation (|12|). Expanding all exponents in (|12|) and performing the 
functional differentiations results in the quantity $ being expressed as a diagram series. We 
see that the causal regularisation is exactly what is necessary to eliminate the "short-circuited" 
diagrams (containing D-r(r,f,0) which is undefined), giving rise to infinities. These diagrams 
emerge when the differential operator A is applied to a single 7i, whereas their absence is 
required by the normal form of the interaction (no contractions within the same interaction 
Hamiltonian) . Note that this regularisation makes sense of the replacement £<^(0) — > £tpQ we 
used above. 

The causal regularisation of Dt is also what leads to the interpretation of (|15| ) as an 
Ito SDE. Since the noise source r)(r,T) in ( |15| ) is singular (delta-correlated white noise), 
the SDE is mathematically speaking not defined. However the corresponding integral equa- 
tion can consistently be interpreted H with a proper definition of the stochastic incre- 
ment, ip(T)r](T)dT = ip(r)dW, where W{r) = J 77(r')dr'. In the Ito stochastic calculus 
ip(T)dW is understood as iP(t)[W(t + dr) — W(t)], while in Stratonovich calculus it is 
|[^(t) + ip(r + dr)][W{T + dr) - W(t)]. Alternatively JL3J, one may replace (]l|) by the 
corresponding integral equation, ip = ip + iDt^, which is defined given Dt is regularised. 
It is then easy to see that the causal regularisation of Dt results in ip(r) being uncorrelated 
with T)(t), which is the characteristic property of the Ito calculus. More detailed arguments 
in favour of this conjecture are given in p3|. 

It is worth noting that the diagram series for $ is structurally identical to that for Bose- 
condensed systems [0, pf , where the propagator is replaced by Dt and the condensate amplitude 
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by the coherent amplitude tpQ. Yet, despite this formal similarity, these series are clearly 
distinct. Here the coherent amplitude ipo applies to all bosonic states, not just the condensate. 
This makes it unnecessary to treat condensed and non-condensed fractions separately as in 
[pi p|. A novel feature is furthermore that our series as a whole and not the propagator is 
averaged over the free distribution -Pf^o)- Consequently we find the propagator Dt equal to 
the vacuum average (|10|), whereas in |2j, y] it is a thermal average. 



Quantum oscillator. - To verify relation ( |l8| ) including the said regularisation/calculus 
conjecture, consider an anharmonic oscillator described by annihilation and creation operators 
if) and ip> with H n = ujo'ip^'ip an d ^int = K /2 ^ 2 ip 2 . The above results apply by simply 
dropping the spatial variable, ip(f,T) — » iI>(t), etc. Then, V'o( T ) = ipo e ~ UJ ° T \ V>o( T ) = ^oe UJ ° T . 
The SDE. which now reads 



— - lu - ir]{r) 



$<j) = 0, (20) 



is readily solved Q, which yields, ^i(O) = ztp , where z = se l ® , ■& — /^dr^r) is the total 

random increment, d 2 — (3k, and s — e@ K ' 2 . This result corresponds to (|2Q) regarded as an Ito 
equation. In Stratonovich calculus (say), one would find s = 1. Using the explicit expression 
for P(ip ) g, we have (r = e 0u °), 

(r - lje-fr-^l^l'd 2 ^ /" e-**/ 2 ""^ __ r u , 2 , 



*(£,0 = / Z — / ^^ exp[-|^ |^(l-z)+^o-r^o]. (21) 



Changing the order of integrations in (Ell) and then taking the Gaussian integral over tpo yields, 



m,a= ^exp(-i^) =^E L «(iei 2 )(-)V" 2/3K/2 , (22) 

r — z \ r — zj r * — ' Vr/ 



n=0 

where the upper bar denotes the integration (averaging) over d and L„ are the Laguerre 
polynomials. 

To verify the choice of the Ito calculus consider the partition function of the anharmonic 
oscillator. Noting that x(0, 0) = (S) = Z/Z , Z = 1/(1 -e^ ) =r/(r-l), and using @, 
we have 



Z= —^ = £^(-/W»-^Y (23) 



ra=0 



This should be compared to the Fock-space expansion, Z = X) n =o ex P i~ P^o 71 ~ f3Kn(n — l)/2} . 
We see that the "Ito-valued" s — e /3K / 2 is exactly what is required to match (p3J) and this 
expression, while the "Stratonovich-valued" s — 1 results in a discrepancy. 

In Fig.l, we plot the characteristic function for ujq — I, /3 — 0.1 and different values of the 
interaction constant. After an initial decay, \ exhibits exponentially growing oscillations. The 
P-function, which is the Fourier transform of x, is hence not an ordinary function and the 
thermal state of the interacting system is non- classical as opposed to the non- interacting one. 
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Fig. 1. - Characteristic function of the oscillator for uio 
anharmonicity. 



1, f3 — 0.1 and different values of 



The anharmonic oscillator shows an interesting behaviour if a negative linear part is added 
to the interaction, H[ nt = — Auj ffi i[j + k^^ 2 ^ 2 . The total linear part of the Hamiltonian, 
(wo — A(j)h — wft, may then be negative. If ip 2 = —uj/k > 0, the ground state of the 
interacting system no longer coincides with that of the harmonic oscillator but with one of 
the Fock states. In the limit of infinite system-size, hq = ujq/k — > oo, the system undergoes a 
second-order Higgs-type phase transition, with ip being the order parameter. This is illustrated 
in Fig. 2 which shows the normalized mean number of quanta v = (h)/no and fluctuations 
as a function of Alo for different system-size parameters. 




Fig. 2. - Higgs-type phase transition in the anharmonic oscillator. Mean number of quanta and 
fluctuations (insert) for wo = 1, P — 1 and different no- 

Above threshold, quantum distribution such as the Wigner function should become concen- 
trated at \ij}\ ~ ip rather than at \ifj\ ~ 0. The above analysis is easily extended to include the 



additional linear interaction and all results apply with s = e 



/3k/2 



3 /3Ao;+/3 K /2^ Thus 



w(^,r) 



J 7T 7TX(0,0) 



1 



r + z 



exp 



-2|V| 



r + z 



(24) 



The averaging here is over •&. Rigorously the last equation in (p3) holds only if 5R(r — z) > and 
hence k < 2u>, but one can get rid of this condition by appropriately modifying the averaging 
(namely, moving it from the circle \z\ = s to \z\ = 1 while preserving all z n , n = 0, 1, • • ■). 
The results of the corresponding numerical evaluation of the Wigner function shown in Fig. 3 
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clearly indicate the phase transition when lu becomes negative. 





Fig. 3a. - Typical shape of the Wigner function 
"above threshold" (/3 = 0.1, k — 1, w = -8). 



Fig. 3b. 

ues of lu. 



The Wigner function for different val- 



In conclusion, we have shown that in the Feynman-path representation of the Matsubara 
Green's function of interacting bosons, the Feynman paths may be constructively characterised 
as solutions to an imaginary time Ito stochastic differential equation. This allows one to 
calculate normally ordered thermal averages of interacting nonrelativistic Bose fields beyond 
the level of perturbation. We have verified this result, including the fact that the equation 
we find is an Ito equation, for the simple prototype model of an anharmonic oscillator. The 
relation between stochastic calculus and regularisation as well as the application of the method 
to ultracold atomic quantum gases will be subject to further investigations. 
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